clc;
clear;

%% main主脚本
tic;
startTime = cputime;
pde = process_pde_MF(1);
process_plot_ex_scalar(pde, pde.p);
vector_h = [1/2 1/4 1/8];
basis_type1 = "P2";
basis_type2 = "P1";
iter_L = 8;

u_Linf = zeros(size(vector_h,2),1);
u_L2 = zeros(size(vector_h,2),1);
u_H1 = zeros(size(vector_h,2),1);
w_Linf = zeros(size(vector_h,2),1);
w_L2 = zeros(size(vector_h,2),1);
w_H1 = zeros(size(vector_h,2),1);
p_Linf = zeros(size(vector_h,2),1);
p_L2 = zeros(size(vector_h,2),1);
p_H1 = zeros(size(vector_h,2),1);

for i = 1:size(vector_h,2)
    h = vector_h(i);

    mesh = process_mesh(pde, h);
    fs1 = process_fespace(mesh, basis_type1);
    fs2 = process_fespace(mesh, basis_type2);

    % process_showmesh(mesh, fs1);
    % process_showmesh(mesh, fs2);

    A11 = (pde.nu+pde.nu_r)*term_grad_U_grad_V(mesh, fs1);
    A13 = term_P_div_V(mesh, fs1, fs2);
    A14 = 2*pde.nu_r*term_rot_W_V(mesh, fs1);

    A31 = (pde.Ca+pde.Cd)*term_grad_W_grad_Psi(mesh, fs1);
    A33 = 4*pde.nu_r*term_W_Psi(mesh, fs1);

    O1 = sparse(3*fs1.Nb,fs2.Nb);
    O2 = sparse(fs2.Nb,fs2.Nb);
    O3 = sparse(fs2.Nb,1);

    A = [A11, -A14, -A13;
        -A14', A31+A33, O1;
        -A13', O1', O2];

    b1 = term_F_V(mesh, fs1, pde.f11, pde.f12, pde.f13);
    b2 = term_F_V(mesh, fs1, pde.f21, pde.f22, pde.f23);
    b = [b1; b2; O3];

    % U_n = sparse(3*fs1.Nb,1);
    U_n = [pde.u1(fs1.node(:,1), fs1.node(:,2), fs1.node(:,3));
        pde.u2(fs1.node(:,1), fs1.node(:,2), fs1.node(:,3));
        pde.u3(fs1.node(:,1), fs1.node(:,2), fs1.node(:,3))];
    for l = 1:iter_L
        A12 = term_Uold_grad_U_V(mesh, fs1, U_n);
        A32 = term_Uold_grad_W_Psi(mesh, fs1, U_n);
        AN = blkdiag(A12, A32, O2);

        Al = A + AN;
        bl = b;

        [Al, bl] = process_boundary(Al, bl, mesh, fs1, fs2, pde);

        SOL = Al\bl;

        U_n = SOL(1:3*fs1.Nb);
    end

    % process_plot_scalar(mesh, fs1, SOL(1:fs1.Nb));
    % process_plot_scalar(mesh, fs1, SOL(fs1.Nb+1:2*fs1.Nb));
    % process_plot_scalar(mesh, fs1, SOL(2*fs1.Nb+1:3*fs1.Nb));
    % process_plot_scalar(mesh, fs1, SOL(3*fs1.Nb+1:4*fs1.Nb));
    % process_plot_scalar(mesh, fs1, SOL(4*fs1.Nb+1:5*fs1.Nb));
    % process_plot_scalar(mesh, fs1, SOL(5*fs1.Nb+1:6*fs1.Nb));
    process_plot_scalar(mesh, fs2, SOL(6*fs1.Nb+1:end));

    [u_Linf(i), w_Linf(i), p_Linf(i)] = process_error(mesh, fs1, fs2, pde, SOL, "Linf");
    [u_L2(i), w_L2(i), p_L2(i)] = process_error(mesh, fs1, fs2, pde, SOL, "L2");
    [u_H1(i), w_H1(i), p_H1(i)] = process_error(mesh, fs1, fs2, pde, SOL, "H1");

    fprintf("已完成\t%d/%d\n", i, size(vector_h,2));
end

process_rate(u_Linf, vector_h);
process_rate(u_L2, vector_h);
process_rate(u_H1, vector_h);
process_rate(w_Linf, vector_h);
process_rate(w_L2, vector_h);
process_rate(w_H1, vector_h);
process_rate(p_Linf, vector_h);
process_rate(p_L2, vector_h);
process_rate(p_H1, vector_h);

TOCtime = toc;
CPUtime = cputime - startTime;
disp(["TOC time(s)", num2str(TOCtime)]);
disp(["CPU time(s)", num2str(CPUtime)]);




%% process_pde_MF
function pde = process_pde_MF(index)
switch index
    case 1
        pde.index = index;
        nu = 1;
        nu_r = 1;
        C0 = 1;
        Ca = 1;
        Cd = 1;
        pde.nu = nu;
        pde.nu_r = nu_r;
        pde.C0 = C0;
        pde.Ca = Ca;
        pde.Cd = Cd;
        pde.u1 = @(x,y,z)y.^4+z.^2;
        pde.u2 = @(x,y,z)z.^4+x.^2;
        pde.u3 = @(x,y,z)x.^4+y.^2;
        pde.w1 = @(x,y,z)sin(y)+z;
        pde.w2 = @(x,y,z)sin(z)+x;
        pde.w3 = @(x,y,z)sin(x)+y;
        pde.p = @(x,y,z)(2*x-1).*(2*y-1).*(2*z-1);
        pde.u1dx = @(x,y,z)0.0;
        pde.u1dy = @(x,y,z)y.^3.*4.0;
        pde.u1dz = @(x,y,z)z.*2.0;
        pde.u2dx = @(x,y,z)x.*2.0;
        pde.u2dy = @(x,y,z)0.0;
        pde.u2dz = @(x,y,z)z.^3.*4.0;
        pde.u3dx = @(x,y,z)x.^3.*4.0;
        pde.u3dy = @(x,y,z)y.*2.0;
        pde.u3dz = @(x,y,z)0.0;
        pde.w1dx = @(x,y,z)0.0;
        pde.w1dy = @(x,y,z)cos(y);
        pde.w1dz = @(x,y,z)1.0;
        pde.w2dx = @(x,y,z)1.0;
        pde.w2dy = @(x,y,z)0.0;
        pde.w2dz = @(x,y,z)cos(z);
        pde.w3dx = @(x,y,z)cos(x);
        pde.w3dy = @(x,y,z)1.0;
        pde.w3dz = @(x,y,z)0.0;
        pde.pdx = @(x,y,z)(y.*2.0-1.0).*(z.*2.0-1.0).*2.0;
        pde.pdy = @(x,y,z)(x.*2.0-1.0).*(z.*2.0-1.0).*2.0;
        pde.pdz = @(x,y,z)(x.*2.0-1.0).*(y.*2.0-1.0).*2.0;
        pde.f11 = @(x,y,z)nu_r.*(cos(z)-1.0).*2.0-(nu+nu_r).*(y.^2.*1.2e+1+2.0)+y.^3.*(x.^2+z.^4).*4.0+(y.*2.0-1.0).*(z.*2.0-1.0).*2.0+z.*(x.^4+y.^2).*2.0;
        pde.f12 = @(x,y,z)nu_r.*(cos(x)-1.0).*2.0-(nu+nu_r).*(z.^2.*1.2e+1+2.0)+z.^3.*(x.^4+y.^2).*4.0+(x.*2.0-1.0).*(z.*2.0-1.0).*2.0+x.*(y.^4+z.^2).*2.0;
        pde.f13 = @(x,y,z)nu_r.*(cos(y)-1.0).*2.0-(nu+nu_r).*(x.^2.*1.2e+1+2.0)+x.^3.*(y.^4+z.^2).*4.0+(x.*2.0-1.0).*(y.*2.0-1.0).*2.0+y.*(x.^2+z.^4).*2.0;
        pde.f21 = @(x,y,z)sin(y).*(Ca+Cd)+cos(y).*(x.^2+z.^4)+nu_r.*(z+sin(y)).*4.0+x.^4+y.^2-nu_r.*(y-z.^3.*2.0).*4.0;
        pde.f22 = @(x,y,z)sin(z).*(Ca+Cd)+cos(z).*(x.^4+y.^2)+nu_r.*(x+sin(z)).*4.0+y.^4+z.^2-nu_r.*(z-x.^3.*2.0).*4.0;
        pde.f23 = @(x,y,z)sin(x).*(Ca+Cd)+cos(x).*(y.^4+z.^2)+nu_r.*(y+sin(x)).*4.0+x.^2+z.^4-nu_r.*(x-y.^3.*2.0).*4.0;
        % pde.domain = [-0.5 0.5 -0.5 0.5 -0.5 0.5];
        pde.domain = [0 1 0 1 0 1];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax,zmin,zmax\n");
disp(pde.domain);
if isfield(pde, "bdary_u")
    fprintf("pde boundary: left,front,right,back,bottom,top (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_u);
end
if isfield(pde, "bdary_B")
    fprintf("pde boundary of B: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_B);
end
if isfield(pde, "bdary_w")
    fprintf("pde boundary of w: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_w);
end
if isfield(pde, "bdary_p")
    fprintf("pde boundary of P: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_p);
end
end

%% term_grad_U_grad_V
function matrix_A = term_grad_U_grad_V(mesh, fespace)
A1 = term_grad_P_grad_Q(mesh, fespace);
matrix_A = blkdiag(A1, A1, A1);
end

%% term_P_div_V
function matrix_A = term_P_div_V(mesh, fespace1, fespace2)
Nb1 = fespace1.Nb;
Nlb1 = fespace1.Nlb;
Nb2 = fespace2.Nb;
Nlb2 = fespace2.Nlb;
rows = zeros(mesh.Ne*Nlb1*Nlb2,1);
cols = rows;
val1 = rows;
val2 = rows;
val3 = rows;
flag = 1;
for n = 1:mesh.Ne
    V = mesh.node(mesh.elem(n,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dy");
    phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dz");
    psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "x");
    for i = 1:Nlb1
        for j = 1:Nlb2
            f1 = psi(:,j) .* phix(:,i);
            f2 = psi(:,j) .* phiy(:,i);
            f3 = psi(:,j) .* phiz(:,i);
            r1 = int_weights_loc * f1;
            r2 = int_weights_loc * f2;
            r3 = int_weights_loc * f3;
            rows(flag) = fespace1.elem(n,i);
            cols(flag) = fespace2.elem(n,j);
            val1(flag) = r1;
            val2(flag) = r2;
            val3(flag) = r3;
            flag = flag + 1;
        end
    end
end
A1 = sparse(rows, cols, val1, Nb1, Nb2);
A2 = sparse(rows, cols, val2, Nb1, Nb2);
A3 = sparse(rows, cols, val3, Nb1, Nb2);
matrix_A = [A1; A2; A3];
end

%% term_rot_W_V
function matrix_A = term_rot_W_V(mesh, fespace)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb*Nlb,1);
cols = rows;
val1 = rows;
val2 = rows;
val3 = rows;
flag = 1;
for n = 1:mesh.Ne
    V = mesh.node(mesh.elem(n,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
    phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "x");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dy");
    phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dz");
    for i = 1:Nlb
        for j = 1:Nlb
            f1 = phix(:,j).*phi(:,i);
            f2 = phiy(:,j).*phi(:,i);
            f3 = phiz(:,j).*phi(:,i);
            r1 = int_weights_loc * f1;
            r2 = int_weights_loc * f2;
            r3 = int_weights_loc * f3;
            rows(flag) = fespace.elem(n,i);
            cols(flag) = fespace.elem(n,j);
            val1(flag) = r1;
            val2(flag) = r2;
            val3(flag) = r3;
            flag = flag + 1;
        end
    end
end
Ax = sparse(rows, cols, val1, Nb, Nb);
Ay = sparse(rows, cols, val2, Nb, Nb);
Az = sparse(rows, cols, val3, Nb, Nb);
O = sparse(Nb,Nb);
matrix_A = [O, -Az, Ay;
    Az, O, -Ax;
    -Ay, Ax, O];
end

%% term_grad_W_grad_Psi
function matrix_A = term_grad_W_grad_Psi(mesh, fespace)
matrix_A = term_grad_U_grad_V(mesh, fespace);
end

%% term_W_Psi
function matrix_A = term_W_Psi(mesh, fespace)
A1 = term_P_Q(mesh, fespace);
matrix_A = blkdiag(A1, A1, A1);
end

%% term_F_V
function vector_b = term_F_V(mesh, fespace, func_f1, func_f2, func_f3)
b1 = term_F_Psi(mesh, fespace, func_f1);
b2 = term_F_Psi(mesh, fespace, func_f2);
b3 = term_F_Psi(mesh, fespace, func_f3);
vector_b = [b1; b2; b3];
end

%% term_Uold_grad_U_V
function matrix_A = term_Uold_grad_U_V(mesh, fespace, U_n)
Nb = fespace.Nb;
Nlb = fespace.Nlb;
rows = zeros(mesh.Ne*Nlb*Nlb,1);
cols = rows;
val = rows;
flag = 1;
for e = 1:mesh.Ne
    V = mesh.node(mesh.elem(e,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
    phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "x");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dy");
    phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace.basis_type, "dz");
    u1 = phi * U_n(fespace.elem(e,:));
    u2 = phi * U_n(fespace.elem(e,:)+Nb);
    u3 = phi * U_n(fespace.elem(e,:)+2*Nb);
    for i = 1:Nlb
        for j = 1:Nlb
            f = u1.*phix(:,j).*phi(:,i) + u2.*phiy(:,j).*phi(:,i) + u3.*phiz(:,j).*phi(:,i);
            r = int_weights_loc * f;
            rows(flag) = fespace.elem(e,i);
            cols(flag) = fespace.elem(e,j);
            val(flag) = r;
            flag = flag + 1;
        end
    end
end
A1 = sparse(rows, cols, val, Nb, Nb);
matrix_A = blkdiag(A1, A1, A1);
end

%% term_Uold_grad_W_Psi
function matrix_A = term_Uold_grad_W_Psi(mesh, fespace, U_n)
matrix_A = term_Uold_grad_U_V(mesh, fespace, U_n);
end

%% process_boundary
function [A, b] = process_boundary(A, b, mesh, fespace1, fespace2, pde)
Nb1 = fespace1.Nb;
% Neumann

% Dirichlet
for k = 1:6
    % u1
    bdn1 = bdary_node(mesh, fespace1, k);
    A(bdn1,:) = 0;
    A(sub2ind(size(A), bdn1, bdn1)) = 1;
    b(bdn1) = pde.u1(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
    % u2
    bdn2 = bdn1 + Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.u2(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
    % u3
    bdn2 = bdn1 + 2*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.u3(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
    % w1
    bdn2 = bdn1 + 3*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.w1(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
    % w2
    bdn2 = bdn1 + 4*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.w2(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
    % w3
    bdn2 = bdn1 + 5*Nb1;
    A(bdn2,:) = 0;
    A(sub2ind(size(A), bdn2, bdn2)) = 1;
    b(bdn2) = pde.w3(fespace1.node(bdn1,1),fespace1.node(bdn1,2),fespace1.node(bdn1,3));
end
% p
A(6*Nb1+1,:) = 0;
A(6*Nb1+1,6*Nb1+1) = 1;
b(6*Nb1+1,1) = pde.p(pde.domain(1), pde.domain(3), pde.domain(5));
% M = term_P_Q(mesh, fespace2);
% matrix_A = blkdiag(sparse(6*Nb1,6*Nb1), -1e-6*M);
% A = A + matrix_A;
end

%% process_error
function [error_u, error_w, error_p] = process_error(mesh, fespace1, fespace2, pde, SOL, error_type)
Nb1 = fespace1.Nb;
switch error_type
    case "Linf"
        err_u1 = 0;
        err_u2 = 0;
        err_u3 = 0;
        err_w1 = 0;
        err_w2 = 0;
        err_w3 = 0;
        err_p = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, ~] = integral_loc(V, "3d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "x");
            u1h = phi * SOL(fespace1.elem(e,:));
            u2h = phi * SOL(fespace1.elem(e,:)+Nb1);
            u3h = phi * SOL(fespace1.elem(e,:)+2*Nb1);
            w1h = phi * SOL(fespace1.elem(e,:)+3*Nb1);
            w2h = phi * SOL(fespace1.elem(e,:)+4*Nb1);
            w3h = phi * SOL(fespace1.elem(e,:)+5*Nb1);
            ph = psi * SOL(fespace2.elem(e,:)+6*Nb1);
            f1 = max(abs(pde.u1(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u1h));
            f2 = max(abs(pde.u2(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u2h));
            f3 = max(abs(pde.u3(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u3h));
            f4 = max(abs(pde.w1(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w1h));
            f5 = max(abs(pde.w2(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w2h));
            f6 = max(abs(pde.w3(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w3h));
            f7 = max(abs(pde.p(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - ph));
            if f1 > err_u1
                err_u1 = f1;
            end
            if f2 > err_u2
                err_u2 = f2;
            end
            if f3 > err_u3
                err_u3 = f3;
            end
            if f4 > err_w1
                err_w1 = f4;
            end
            if f5 > err_w2
                err_w2 = f5;
            end
            if f6 > err_w3
                err_w3 = f6;
            end
            if f7 > err_p
                err_p = f7;
            end
        end
        error_u = max([err_u1, err_u2, err_u3]);
        error_w = max([err_w1, err_w2, err_w3]);
        error_p = err_p;
    case "L2"
        err_u1 = 0;
        err_u2 = 0;
        err_u3 = 0;
        err_w1 = 0;
        err_w2 = 0;
        err_w3 = 0;
        err_p = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "x");
            psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "x");
            u1h = phi * SOL(fespace1.elem(e,:));
            u2h = phi * SOL(fespace1.elem(e,:)+Nb1);
            u3h = phi * SOL(fespace1.elem(e,:)+2*Nb1);
            w1h = phi * SOL(fespace1.elem(e,:)+3*Nb1);
            w2h = phi * SOL(fespace1.elem(e,:)+4*Nb1);
            w3h = phi * SOL(fespace1.elem(e,:)+5*Nb1);
            ph = psi * SOL(fespace2.elem(e,:)+6*Nb1);
            f1 = (pde.u1(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u1h).^2;
            f2 = (pde.u2(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u2h).^2;
            f3 = (pde.u3(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u3h).^2;
            f4 = (pde.w1(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w1h).^2;
            f5 = (pde.w2(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w2h).^2;
            f6 = (pde.w3(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w3h).^2;
            f7 = (pde.p(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - ph).^2;
            r1 = int_weights_loc * f1;
            r2 = int_weights_loc * f2;
            r3 = int_weights_loc * f3;
            r4 = int_weights_loc * f4;
            r5 = int_weights_loc * f5;
            r6 = int_weights_loc * f6;
            r7 = int_weights_loc * f7;
            err_u1 = err_u1 + r1;
            err_u2 = err_u2 + r2;
            err_u3 = err_u3 + r3;
            err_w1 = err_w1 + r4;
            err_w2 = err_w2 + r5;
            err_w3 = err_w3 + r6;
            err_p = err_p + r7;
        end
        error_u = sqrt(err_u1 + err_u2 + err_u3);
        error_w = sqrt(err_w1 + err_w2 + err_w3);
        error_p = sqrt(err_p);
    case "H1"
        err_u1x = 0;
        err_u1y = 0;
        err_u1z = 0;
        err_u2x = 0;
        err_u2y = 0;
        err_u2z = 0;
        err_u3x = 0;
        err_u3y = 0;
        err_u3z = 0;
        err_w1x = 0;
        err_w1y = 0;
        err_w1z = 0;
        err_w2x = 0;
        err_w2y = 0;
        err_w2z = 0;
        err_w3x = 0;
        err_w3y = 0;
        err_w3z = 0;
        err_px = 0;
        err_py = 0;
        err_pz = 0;
        for e = 1:mesh.Ne
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(V, "3d");
            phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dx");
            phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dy");
            phiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace1.basis_type, "dz");
            psix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "dx");
            psiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "dy");
            psiz = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), int_points_loc(:,3), fespace2.basis_type, "dz");
            u1hx = phix * SOL(fespace1.elem(e,:));
            u1hy = phiy * SOL(fespace1.elem(e,:));
            u1hz = phiz * SOL(fespace1.elem(e,:));
            u2hx = phix * SOL(fespace1.elem(e,:)+Nb1);
            u2hy = phiy * SOL(fespace1.elem(e,:)+Nb1);
            u2hz = phiz * SOL(fespace1.elem(e,:)+Nb1);
            u3hx = phix * SOL(fespace1.elem(e,:)+2*Nb1);
            u3hy = phiy * SOL(fespace1.elem(e,:)+2*Nb1);
            u3hz = phiz * SOL(fespace1.elem(e,:)+2*Nb1);
            w1hx = phix * SOL(fespace1.elem(e,:)+3*Nb1);
            w1hy = phiy * SOL(fespace1.elem(e,:)+3*Nb1);
            w1hz = phiz * SOL(fespace1.elem(e,:)+3*Nb1);
            w2hx = phix * SOL(fespace1.elem(e,:)+4*Nb1);
            w2hy = phiy * SOL(fespace1.elem(e,:)+4*Nb1);
            w2hz = phiz * SOL(fespace1.elem(e,:)+4*Nb1);
            w3hx = phix * SOL(fespace1.elem(e,:)+5*Nb1);
            w3hy = phiy * SOL(fespace1.elem(e,:)+5*Nb1);
            w3hz = phiz * SOL(fespace1.elem(e,:)+5*Nb1);
            phx = psix * SOL(fespace2.elem(e,:)+6*Nb1);
            phy = psiy * SOL(fespace2.elem(e,:)+6*Nb1);
            phz = psiz * SOL(fespace2.elem(e,:)+6*Nb1);
            f1x = (pde.u1dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u1hx).^2;
            f1y = (pde.u1dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u1hy).^2;
            f1z = (pde.u1dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u1hz).^2;
            f2x = (pde.u2dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u2hx).^2;
            f2y = (pde.u2dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u2hy).^2;
            f2z = (pde.u2dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u2hz).^2;
            f3x = (pde.u3dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u3hx).^2;
            f3y = (pde.u3dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u3hy).^2;
            f3z = (pde.u3dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - u3hz).^2;
            f4x = (pde.w1dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w1hx).^2;
            f4y = (pde.w1dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w1hy).^2;
            f4z = (pde.w1dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w1hz).^2;
            f5x = (pde.w2dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w2hx).^2;
            f5y = (pde.w2dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w2hy).^2;
            f5z = (pde.w2dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w2hz).^2;
            f6x = (pde.w3dx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w3hx).^2;
            f6y = (pde.w3dy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w3hy).^2;
            f6z = (pde.w3dz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - w3hz).^2;
            f7x = (pde.pdx(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - phx).^2;
            f7y = (pde.pdy(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - phy).^2;
            f7z = (pde.pdz(int_points_loc(:,1),int_points_loc(:,2),int_points_loc(:,3)) - phz).^2;
            r1x = int_weights_loc * f1x;
            r1y = int_weights_loc * f1y;
            r1z = int_weights_loc * f1z;
            r2x = int_weights_loc * f2x;
            r2y = int_weights_loc * f2y;
            r2z = int_weights_loc * f2z;
            r3x = int_weights_loc * f3x;
            r3y = int_weights_loc * f3y;
            r3z = int_weights_loc * f3z;
            r4x = int_weights_loc * f4x;
            r4y = int_weights_loc * f4y;
            r4z = int_weights_loc * f4z;
            r5x = int_weights_loc * f5x;
            r5y = int_weights_loc * f5y;
            r5z = int_weights_loc * f5z;
            r6x = int_weights_loc * f6x;
            r6y = int_weights_loc * f6y;
            r6z = int_weights_loc * f6z;
            r7x = int_weights_loc * f7x;
            r7y = int_weights_loc * f7y;
            r7z = int_weights_loc * f7z;
            err_u1x = err_u1x + r1x;
            err_u1y = err_u1y + r1y;
            err_u1z = err_u1z + r1z;
            err_u2x = err_u2x + r2x;
            err_u2y = err_u2y + r2y;
            err_u2z = err_u2z + r2z;
            err_u3x = err_u3x + r3x;
            err_u3y = err_u3y + r3y;
            err_u3z = err_u3z + r3z;
            err_w1x = err_w1x + r4x;
            err_w1y = err_w1y + r4y;
            err_w1z = err_w1z + r4z;
            err_w2x = err_w2x + r5x;
            err_w2y = err_w2y + r5y;
            err_w2z = err_w2z + r5z;
            err_w3x = err_w3x + r6x;
            err_w3y = err_w3y + r6y;
            err_w3z = err_w3z + r6z;
            err_px = err_px + r7x;
            err_py = err_py + r7y;
            err_pz = err_pz + r7z;
        end
        error_u = sqrt(err_u1x + err_u1y + err_u1z + err_u2x + err_u2y + err_u2z + err_u3x + err_u3y + err_u3z);
        error_w = sqrt(err_w1x + err_w1y + err_w1z + err_w2x + err_w2y + err_w2z + err_w3x + err_w3y + err_w3z);
        error_p = sqrt(err_px + err_py + err_pz);
    otherwise
        error("Invalid error type.");
end
end